/* Bessel I_1(x) function in double precision */

#include <math.h>

double dbesi1(double x)
{
    int k;
    double w, t, y;
    static double a[60] = {
        1.2787464404046789181e-10, 3.5705860060088241077e-9, 
        9.961153761934733504e-8, 2.2395070088633043177e-6, 
        4.0312466928887462346e-5, 5.6437387840203722356e-4, 
        0.0059259259312934746096, 0.04444444444349900887, 
        0.22222222222232042719, 0.66666666666666139867, 
        1.0000000000000001106, 0.49999999999999999962, 
        1.7281952384448634449e-10, 3.064720455997639013e-9, 
        1.0237662138842827028e-7, 2.2299494417341498163e-6, 
        4.0335364374929326943e-5, 5.6433440269141349899e-4, 
        0.0059259754885893798654, 0.04444439941088039787, 
        0.2222222511283502673, 0.66666665422146063244, 
        1.0000000032274936821, 0.49999999961866867205, 
        2.3216048939948030996e-10, 1.7443372702334489579e-9, 
        1.1596478963485415499e-7, 2.1446755518623035147e-6, 
        4.0697440347437076195e-5, 5.6324394900433192204e-4, 
        0.0059283484996093060678, 0.044440673899150997921, 
        0.2222263801685265786, 0.66666358151576732094, 
        1.0000013834029985337, 0.49999971643129650249, 
        3.1013758938255172562e-10, -8.4813676145611694984e-10, 
        1.5544980187411802596e-7, 1.7811109378708045726e-6, 
        4.2945322199060856985e-5, 5.5344850176852353639e-4, 
        0.0059590327716950614802, 0.044371611097707060659, 
        0.22233578241986401111, 0.6665474730046331531, 
        1.0000756505206705927, 0.49997803664415994554, 
        4.1214758313965020365e-10, -5.361331773534742944e-9, 
        2.4661360807517345161e-7, 6.7144593918926723203e-7, 
        5.1988027944945587571e-5, 5.0165568586065803067e-4, 
        0.0061717530047005289953, 0.043745229577317251404, 
        0.22363147971477747996, 0.6647546913111766024, 
        1.0015686689447547657, 0.49941120439785391891
    };
    static double b[70] = {
        6.6324787943143095845e-8, 4.5125928898466638619e-7, 
        6.7937793134877246623e-6, 7.4580507871505926302e-5, 
        7.6866382927334005919e-4, 0.0071185174803491859307, 
        0.058721838073486424416, 0.42473949281714196041, 
        2.6396965606282079123, 13.710008536637016903, 
        57.158647688180932003, 179.46182892089389037, 
        377.57997362398478619, 399.87313678256009819, 
        1.7652713206027939711e-7, 1.1988179244834708057e-6, 
        1.8037851545747139231e-5, 1.9775785516370314656e-4, 
        0.0020354870702829387283, 0.0188221641910322536, 
        0.15500485219010424263, 1.119010001056057321, 
        6.9391565185406617552, 35.948170579648649345, 
        149.41909525103032616, 467.42979492780642582, 
        979.04227423171290408, 1030.9147225169564443, 
        4.7022299276154507603e-7, 3.1878571710170115972e-6, 
        4.7940153875711448496e-5, 5.2496623508411440227e-4, 
        0.0053968661134780824779, 0.049837081920693776234, 
        0.40979593830387765545, 2.9533186922862948404, 
        18.278176130722516369, 94.47649715018912107, 
        391.66075612645333624, 1221.4182034643210345, 
        2548.6177980961291004, 2670.9883037012546541, 
        1.2535083724002034147e-6, 8.484587142065570825e-6, 
        1.2753227372734042108e-4, 0.0013950105363562648921, 
        0.014325473993765291906, 0.13212452778932829125, 
        1.0849287786885151432, 7.8068089156260172673, 
        48.232254570679165833, 248.80659424902394371, 
        1029.0736929484210803, 3200.5629438795801652, 
        6656.7749162019607914, 6948.8586598121632302, 
        3.3439394490599745013e-6, 2.2600596902211837757e-5, 
        3.3955927589987356838e-4, 0.0037105306061050972474, 
        0.038065263634919156421, 0.35068223415665236079, 
        2.8760027832105027316, 20.665999500843274339, 
        127.47939148516390205, 656.43636874254000885, 
        2709.524283793247992, 8407.1174233600734871, 
        17437.146284159740233, 18141.3487816388316
    };
    static double c[45] = {
        -2.8849790431465382128e-15, -3.5125350943844774657e-14, 
        -7.485086701370741975e-13, -1.8383904048277485153e-11, 
        -5.7303556446977223342e-10, -2.4449502737311496525e-8, 
        -1.6765373351766929724e-6, -3.2189516835265773471e-4, 
        0.051503226936425277377, 
        -5.8674306822281631119e-15, -9.4884898451194085565e-15, 
        -8.503386513660036434e-13, -1.8142997866945285736e-11, 
        -5.7340238386338193949e-10, -2.4449138101742183665e-8, 
        -1.6765375646678855842e-6, -3.2189516826945356325e-4, 
        0.051503226936412017608, 
        -1.4723362506764340882e-14, 1.3945147385179042899e-13, 
        -1.9618041857586930923e-12, -1.3343606394065121821e-11, 
        -5.8649674606973244159e-10, -2.4426060539669553778e-8, 
        -1.6765631828366988006e-6, -3.2189515191449587253e-4, 
        0.051503226931820146445, 
        -5.8203519372580372987e-14, 1.2266326995309845825e-12, 
        -1.3921625844526453237e-11, 6.2228025878281625469e-11, 
        -8.8636681342142794023e-10, -2.3661241616744818608e-8, 
        -1.6777870960740520557e-6, -3.2189402882677074318e-4, 
        0.051503226479551959376, 
        -4.5801527369223291722e-14, 6.7998819697143727209e-13, 
        -4.1624857909290468421e-12, -3.2849009406112440998e-11, 
        -3.247827569043111827e-10, -2.5739209934053714983e-8, 
        -1.6730566573215739195e-6, -3.2190010909008684076e-4, 
        0.05150322986693207715
    };

    w = fabs(x);
    if (w < 8.5) {
        t = w * w * 0.0625;
        k = 12 * ((int) t);
        y = (((((((((((a[k] * t + a[k + 1]) * t + 
            a[k + 2]) * t + a[k + 3]) * t + a[k + 4]) * t + 
            a[k + 5]) * t + a[k + 6]) * t + a[k + 7]) * t + 
            a[k + 8]) * t + a[k + 9]) * t + a[k + 10]) * t + 
            a[k + 11]) * w;
    } else if (w < 12.5) {
        k = (int) w;
        t = w - k;
        k = 14 * (k - 8);
        y = ((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * t + b[k + 13];
    } else {
        t = 60 / w;
        k = 9 * ((int) t);
        y = ((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * sqrt(t) * exp(w);
    }
    return x < 0 ? -y : y;
}

